## this script analyzes the individual-level differences ##
## in voter reactions to compromise ##

	library(foreign)
	setwd('replication')

## create model objects ##

	model1 <- list()
	model2 <- list()
	
## let's go ##

	for(i in 1:10){
		
		data <- read.dta(paste('mixBetas', i, '.dta', sep = ''))
		data$pec <- 0
		data$pec[data$survey == 'gles02.05'] <- 1
		data$pec[data$survey == 'dpes86.89'] <- 1
		data$pec[data$survey == 'dpes89.94'] <- 1
		data$pec[data$survey == 'nes01.05'] <- 1
		data$pec[data$survey == 'ses91.94'] <- 1

		bbt <- na.omit(read.dta(paste('fullBeta', i, '.dta', sep = '')))

		bbt <- colMeans(bbt)

		data$govCompSubNew <- c()
		data$govCompSubNew[data$survey == 'dpes82.86'] <- sd(data$govCompSub[data$survey == 'dpes82.86' & data$gov == 1])
		data$govCompSubNew[data$survey == 'dpes86.89'] <- sd(data$govCompSub[data$survey == 'dpes86.89' & data$gov == 1])
		data$govCompSubNew[data$survey == 'dpes89.94'] <- sd(data$govCompSub[data$survey == 'dpes89.94' & data$gov == 1])
		data$govCompSubNew[data$survey == 'nes01.05'] <- sd(data$govCompSub[data$survey == 'nes01.05' & data$gov == 1])
		data$govCompSubNew[data$survey == 'nzes05.08'] <- sd(data$govCompSub[data$survey == 'nzes05.08' & data$gov == 1])
		data$govCompSubNew[data$survey == 'dnes01.05'] <- sd(data$govCompSub[data$survey == 'dnes01.05' & data$gov == 1])
		data$govCompSubNew[data$survey == 'ses91.94'] <- sd(data$govCompSub[data$survey == 'ses91.94' & data$gov == 1])
		data$govCompSubNew[data$survey == 'gles02.05'] <- sd(data$govCompSub[data$survey == 'gles02.05' & data$gov == 1])
		data$govCompSubNew[data$survey == 'gles05.09'] <- sd(data$govCompSub[data$survey == 'gles05.09' & data$gov == 1])
		data$govCompSubOld <- (data$govCompSub - data$govCompSubNew) * data$gov
		data$govCompSubNew <- (data$govCompSub + data$govCompSubNew) * data$gov

		data$holdL <- exp(data$randGov * data$gov + 
			data$randComp * data$govCompSubOld +
			data$randLr * data$lr +
			data$ranLrgGov * (data$lr * data$gov) + 
			data$randEcon * data$economyGov + 
			data$randDisc * data$govCompSubOld * data$lr * data$gov + 
			data$randDiscEcon * data$govCompSubOld * data$economyGov)

		data$changeL <- exp(data$randGov * data$gov + 
			data$randComp * data$govCompSubOld +
			data$randLr * data$lr +
			data$ranLrgGov * (data$lr * data$gov) + 
			data$randEcon * data$economyGov + 
			data$randDisc * data$govCompSubNew * data$lr * data$gov + 
			data$randDiscEcon * data$govCompSubNew * data$economyGov)
						
		respExtreme <- c()
		polIntSelf1 <- c()
		partisan <- c()
		duration <- c()
		priorGovSupporter <- c()
		survey <- c()
		randGov <- c()
		randComp <- c()
		randLr <- c()
		ranLrgGov <- c()
		randEcon <- c()
		randDisc <- c()
		randDiscEcon <- c()
		probHold <- c()
		probChange <- c()
		pec <- c()		
		survey <- c()		
		change <- c()
		
		group <- sort(unique(data$group))

		for(k in group){
			respExtreme[k] <- mean(data$respExtreme[data$group == k], na.rm = TRUE)
			polIntSelf1[k] <- mean(data$polIntSelf1[data$group == k], na.rm = TRUE)
			partisan[k] <- max(data$partisanStrength1[data$group == k], na.rm = TRUE)
			priorGovSupporter[k] <- mean(data$priorGovSupporter[data$group == k], na.rm = TRUE)
			pec[k] <- data$pec[data$group == k][1]
			survey[k] <- as.character(data$survey[data$group == k][1])

			probHold[k] <- sum(data$holdL[data$gov == 1 & data$group == k])/sum(data$holdL[data$group == k])
			probChange[k] <- sum(data$changeL[data$gov == 1 & data$group == k])/sum(data$changeL[data$group == k])
			change[k] <- probChange[k] - probHold[k]
		}
		
		respExtreme[survey == 'gles02.05'] <- scale(respExtreme[survey == 'gles02.05'])
		respExtreme[survey == 'gles05.09'] <- scale(respExtreme[survey == 'gles05.09'])
		respExtreme[survey == 'dpes82.86'] <- scale(respExtreme[survey == 'dpes82.86'])
		respExtreme[survey == 'dpes86.89'] <- scale(respExtreme[survey == 'dpes86.89'])
		respExtreme[survey == 'dpes89.94'] <- scale(respExtreme[survey == 'dpes89.94'])
		respExtreme[survey == 'nes01.05'] <- scale(respExtreme[survey == 'nes01.05'])
		respExtreme[survey == 'nzes05.08'] <- scale(respExtreme[survey == 'nzes05.08'])
		respExtreme[survey == 'ses91.94'] <- scale(respExtreme[survey == 'ses91.94'])
		respExtreme[survey == 'dnes01.05'] <- scale(respExtreme[survey == 'dnes01.05'])
		
		model1[[i]] <- lm(change ~ respExtreme + polIntSelf1 + partisan + priorGovSupporter + survey)
		summary(model1[[i]])
		model2[[i]] <- lm(change ~ respExtreme + polIntSelf1 + partisan + priorGovSupporter + pec + survey)
		summary(model2[[i]])

		cat(paste('iteration', i, 'complete...\n'))
	
	}

## prediction ##

	zero <- data.frame(respExtreme, polIntSelf1, partisan, priorGovSupporter, survey)
	zero$priorGovSupporter <- 0
	one <- data.frame(respExtreme, polIntSelf1, partisan, priorGovSupporter, survey)
	one$priorGovSupporter <- 1
	
	zero$probs <- predict.lm(model1[[i]], newdata = zero)
	one$probs <- predict.lm(model1[[i]], newdata = one)
	summary(one$probs/zero$probs)
	
	summary(respExtreme)
	sd(respExtreme)
	summary(polIntSelf1)
	sd(polIntSelf1)
	summary(partisan)
	sd(partisan)
	summary(priorGovSupporter)
	sd(priorGovSupporter)
	summary(pec)
	sd(pec)
	gles02.05 <- as.numeric(survey == 'gles02.05')
	gles05.09 <- as.numeric(survey == 'gles05.09')
	dpes82.86 <- as.numeric(survey == 'dpes82.86')
	dpes86.89 <- as.numeric(survey == 'dpes86.89')
	dpes89.94 <- as.numeric(survey == 'dpes89.94')
	nes01.05 <- as.numeric(survey == 'nes01.05')
	nzes05.08 <- as.numeric(survey == 'nzes05.08')
	ses91.94 <- as.numeric(survey == 'ses91.94')
	summary(gles02.05[survey != 'dnes01.05'])
	sd(gles02.05[survey != 'dnes01.05'], na.rm = TRUE)
	summary(gles05.09[survey != 'dnes01.05'])
	sd(gles05.09[survey != 'dnes01.05'], na.rm = TRUE)
	summary(dpes82.86[survey != 'dnes01.05'])
	sd(dpes82.86[survey != 'dnes01.05'], na.rm = TRUE)
	summary(dpes86.89[survey != 'dnes01.05'])
	sd(dpes86.89[survey != 'dnes01.05'], na.rm = TRUE)
	summary(dpes89.94[survey != 'dnes01.05'])
	sd(dpes89.94[survey != 'dnes01.05'], na.rm = TRUE)
	summary(nes01.05[survey != 'dnes01.05'])
	sd(nes01.05[survey != 'dnes01.05'], na.rm = TRUE)
	summary(nzes05.08[survey != 'dnes01.05'])
	sd(nzes05.08[survey != 'dnes01.05'], na.rm = TRUE)
	summary(ses91.94[survey != 'dnes01.05'])
	sd(ses91.94[survey != 'dnes01.05'], na.rm = TRUE)
	
	
## now build the propers SEs, aggregate the fits, and produce the tables ##
## model 1, with no PEC indicator ##

	names <- c('Intercept', 'Respondent Extremity', 'Political Interest', 'Partisanship', 'Prior Cabinet Supporter', 
		'DPES 82-86', 'DPES 86-89', 'DPES 89-94', 'GLES 02-05', 'GLES 05-09', 'NES 01-05', 'NZES 05-08', 'SES 91-94')

	ll <- c()
	coef <- matrix(NA, ncol = 10, nrow = 13)
	se <- matrix(NA, ncol = 10, nrow = 13)
	n <- c()
	r <- c()
	for(i in 1:10){
		ll[i] <- logLik(model1[[i]])
		coef[, i] <- coef(model1[[i]])
		se[, i] <- sqrt(diag(vcov(model1[[i]])))
		n[i] <- as.numeric(summary(model1[[i]][3])[1])
		r[i] <- as.numeric(summary(model1[[i]])[8])
	}
	
## table ##

	par <- c()
	par <- c()
	parVariance <- c()
	meanSeSq <- c()
	seVariance <- c()
	rubin <- c()
	m <- 10
	
	for(i in 1:13){
		par[i] <- mean(coef[i, ])
		parVariance[i] <- sum((par[i] - coef[i, ])^2) / (m - 1)
		meanSeSq[i] <- mean(se[i, ]^2)
		rubin[i] <- sqrt(meanSeSq[i] + (parVariance[i] * (1 + (1/m))))
	}
	
## table for loss model 1 ##

	cbind(names, round(par, 3), round(rubin, 3))
	mean(ll)
	mean(r)
	mean(n)

## now the model with the PECs ##
## ******NOTE****** the Swedish and Kiwi survey fixed effects are collinear and the Swedish indicator is dropped as a result ##

	names <- c('Intercept', 'Respondent Extremity', 'Political Interest', 'Partisanship', 'Prior Cabinet Supporter', 'PEC', 
		'DPES 82-86', 'DPES 86-89', 'DPES 89-94', 'GLES 02-05', 'GLES 05-09', 'NES 01-05', 'NZES 05-08')
	ll <- c()
	coef <- matrix(NA, ncol = 10, nrow = 14)
	se <- matrix(NA, ncol = 10, nrow = 13)
	n <- c()
	r <- c()
	for(i in 1:10){
		ll[i] <- logLik(model2[[i]])
		coef[, i] <- coef(model2[[i]])
		se[, i] <- sqrt(diag(vcov(model2[[i]])))
		n[i] <- as.numeric(summary(model2[[i]][3])[1])
		r[i] <- as.numeric(summary(model2[[i]])[8])
	}
	
	coef <- coef[c(1:13), ]

## table ##

	par <- c()
	par <- c()
	parVariance <- c()
	meanSeSq <- c()
	seVariance <- c()
	rubin <- c()
	m <- 10
	
	for(i in 1:13){
		par[i] <- mean(coef[i, ])
		parVariance[i] <- sum((par[i] - coef[i, ])^2) / (m - 1)
		meanSeSq[i] <- mean(se[i, ]^2)
		rubin[i] <- sqrt(meanSeSq[i] + (parVariance[i] * (1 + (1/m))))
	}
	
## table for loss model 1 ##

	cbind(names, round(par, 3), round(rubin, 3))
	mean(ll)
	mean(r)
	mean(n)